rm(list = ls())

coeftest.cluster <- function(data, fm, cluster1=NULL, cluster2=NULL, ret="test") {
  
  library(sandwich)
  library(lmtest)
  
  data <- as.data.frame(data)
  
  # Return White (1980) standard errors if no cluster
  # variable is provided
  if (is.null(cluster1)) {
    if (ret=="cov") {
      return(vcovHC(fm, type = "HC0"))
    } else {
      return(coeftest(fm, vcov = vcovHC(fm, type = "HC0")))
    }
  }
  
  # Calculation shared by covariance estimates
  est.fun <- estfun(fm)
  
  # Need to identify observations used in the regression (i.e.,
  # non-missing) values, as the cluster vectors come from the full
  # data set and may not be in the regression model.
  inc.obs <- !is.na(est.fun[,1])
  est.fun <- est.fun[inc.obs,]
  
  # Shared data for degrees-of-freedom corrections
  N  <- dim(fm$model)[1]
  NROW <- NROW(est.fun)
  K  <- fm$rank
  
  # Calculate the sandwich covariance estimate
  cov <- function(cluster) {
    cluster <- factor(cluster, exclude=NULL)
    
    # Calculate the "meat" of the sandwich estimators
    u <- apply(est.fun, 2, function(x) tapply(x, cluster, sum))
    meat <- crossprod(u)/N
    
    # Calculations for degrees-of-freedom corrections, followed
    # by calculation of the variance-covariance estimate.
    # NOTE: NROW/N is a kluge to address the fact that sandwich
    # uses the wrong number of rows (includes rows omitted from
    # the regression).
    M <- length(levels(cluster))
    dfc <- M/(M-1) * (N-1)/(N-K)
    
    #print (sandwich(fm, meat=meat))
    return(dfc * NROW/N * sandwich(fm, meat=meat))
  }
  
  # Calculate the covariance matrix estimate for the first cluster.
  cluster1 <- data[inc.obs, cluster1]
  cov1  <- cov(cluster1)
  # print(cov1)
  
  if (is.null(cluster2)) {
    # If only one cluster supplied, return single cluster
    # results
    if (ret=="cov") {
      return(cov1)
    } else {
      return(coeftest(fm, cov1))
    }
  } else {
    # Otherwise do the calculations for the second cluster
    # and the "intersection" cluster.
    cluster2 <- data[inc.obs, cluster2]
    cluster12 <- paste(cluster1, cluster2, sep="")
    
    # Calculate the covariance matrices for cluster2, the "intersection"
    # cluster, then then put all the pieces together.
    cov2   <- cov(cluster2)
    cov12  <- cov(cluster12)
    covMCL <- (cov1 + cov2 - cov12)
    
    # Return the output of coeftest using two-way cluster-robust
    # standard errors.
    # print(ret)
    if (ret=="cov") {
      return(covMCL)
    } else {
      return(coeftest(fm, covMCL))
    }
  }
}


load("dataBJPOLS.RData")

inst.1 <- "judgeiv_hd"
endo.1 <- "pti"
outc.1 <- "vote_post"

time.controls <- "as.factor(court_time1) + as.factor(court_time2) + as.factor(court_dow) + as.factor(severity)"
demo.controls <- "age + I(age^2) +  as.factor(race4) + female + vote_pre + as.factor(noteli) + regis_before"
case.controls <- "as.factor(any_drug) +  as.factor(any_weapon) +  as.factor(any_prop) + as.factor(any_prior_case)"

form.1.fs <- formula(paste(endo.1, "~", inst.1, "+", time.controls))
form.2.fs <- formula(paste(endo.1, "~", inst.1, "+", time.controls, "+", demo.controls))
form.3.fs <- formula(paste(endo.1, "~", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))

form.1 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "|", inst.1, "+", time.controls))
form.2 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "|", inst.1, "+", time.controls, "+", demo.controls))
form.3 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "+", case.controls, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))

m1d1 <- lm(form.1.fs, data = last.cases)
m1d2 <- lm(form.2.fs, data = last.cases)
m1d3 <- lm(form.3.fs, data = last.cases)

set.seed('123')
results.se <- list()
for(i in 1:500) {
  
  out <- last.cases[, this_id_fa[sample.int(.N, .N, TRUE)], by = c("judge_cat", "court_time1")]
  last.cases.boot <- last.cases[last.cases$this_id_fa %in% out$V1, ]
  
  m1d1b <- lm(form.1.fs, data = last.cases.boot)
  m1d2b <- lm(form.2.fs, data = last.cases.boot)
  m1d3b <- lm(form.3.fs, data = last.cases.boot)
 
  results.se[[i]] <- c(coefficients(m1d1b)["judgeiv_hd"],
                       coefficients(m1d2b)["judgeiv_hd"],
                       coefficients(m1d3b)["judgeiv_hd"]) 
}

results.se.f <- apply(do.call('rbind', results.se), 2, sd)

m1a1 <- ivreg(form.1, data = last.cases)
m1a2 <- ivreg(form.2, data = last.cases)
m1a3 <- ivreg(form.3, data = last.cases)

m1a1_d <- summary(m1a1, vcov = sandwich, diagnostic = T)
m1a2_d <- summary(m1a2, vcov = sandwich, diagnostic = T)
m1a3_d <- summary(m1a3, vcov = sandwich, diagnostic = T)

m1a1_d1 <- m1a1_d$diagnostics[grep("Weak", rownames(m1a1_d$diagnostics)), 3]
m1a2_d1 <- m1a2_d$diagnostics[grep("Weak", rownames(m1a2_d$diagnostics)), 3]
m1a3_d1 <- m1a3_d$diagnostics[grep("Weak", rownames(m1a3_d$diagnostics)), 3]

c1 <- round(coeftest.cluster(last.cases, m1d1, cluster1 = "judge_cat"), 2)
c1 <- c1[grep("judgeiv_hd", rownames(c1)), ]

c2 <- round(coeftest.cluster(last.cases, m1d2, cluster1 = "judge_cat"), 2)
c2 <- c2[grep("judgeiv_hd", rownames(c2)), ]

c3 <- round(coeftest.cluster(last.cases, m1d3, cluster1 = "judge_cat"), 2)
c3 <- c3[grep("judgeiv_hd", rownames(c3)), ]

p1 <- data.table(cbind(c1["Estimate"], c2["Estimate"], c3["Estimate"]))
p1$V4 <- 1
p1$V5 <- 1:nrow(p1)
p1$V6 <- "Estimate"

p2 <- data.table(t(as.matrix(results.se.f)))
colnames(p2) <- c("V1", "V2", "V3")
p2$V4 <- 2
p2$V5 <- 1:nrow(p2)
p2$V6 <- "Std. Error"

p3 <- data.table(round(cbind(m1a1_d1, m1a2_d1, m1a3_d1), 2))
colnames(p3) <- c("V1", "V2", "V3")
p3$V4 <- 3
p3$V5 <- 1:nrow(p3)
p3$V6 <- "F-stat"

table <- rbind(p1, p2, p3)
table <- table[order(V5, V4), ]
table$V4 <- table$V5 <- NULL
p4 <- data.table(m1a1$nobs, m1a2$nobs, m1a3$nobs)    
p4$V6 <- "N"

table <- rbind(table, p4)
table <- table[, c("V6", "V1", "V2", "V3")]
colnames(table) <- c("Variable", "Model 1", "Model 2", "Model 3")
table

time.controls <- "as.factor(court_time1) + as.factor(court_time2) + as.factor(court_dow) + as.factor(severity)"
demo.controls <- "age + I(age^2) +  as.factor(race4) + female + vote_pre + as.factor(noteli) + regis_before"
case.controls <- "as.factor(any_drug) +  as.factor(any_weapon) +  as.factor(any_prop) + as.factor(any_prior_case)"

demo.controls_race <- "age + I(age^2)  + female + vote_pre + as.factor(noteli) + regis_before"
demo.controls_gender <- "age + I(age^2) +  as.factor(race4) + vote_pre + as.factor(noteli) + regis_before"

case.controls_drug <- "as.factor(any_weapon) +  as.factor(any_prop) + as.factor(any_prior_case)"
case.controls_prop <- "as.factor(any_weapon) + as.factor(any_drug)  + as.factor(any_prior_case)"
case.controls_weapon <- "as.factor(any_drug) +  as.factor(any_prop) + as.factor(any_prior_case)"
time.controls_violent <- "as.factor(court_time1) + as.factor(court_time2) + as.factor(court_dow)"
case.controls_prior <- "as.factor(any_drug) +  as.factor(any_weapon) +  as.factor(any_prop)"

form.1.fs <- formula(paste(endo.1, "~", inst.1, "+", time.controls))
form.2.fs <- formula(paste(endo.1, "~", inst.1, "+", time.controls, "+", demo.controls))
form.3.fs <- formula(paste(endo.1, "~", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))

form.1 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "|", inst.1, "+", time.controls))
form.2 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "|", inst.1, "+", time.controls, "+", demo.controls))
form.3 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "+", case.controls, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))

form.2.fs_race <- formula(paste(endo.1, "~", inst.1, "+", time.controls, "+", demo.controls_race))
form.3.fs_race <- formula(paste(endo.1, "~", inst.1, "+", time.controls, "+", demo.controls_race, "+", case.controls))

form.2_race <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls_race, "|", inst.1, "+", time.controls, "+", demo.controls_race))
form.3_race <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls_race, "+", case.controls, "|", inst.1, "+", time.controls, "+", demo.controls_race, "+", case.controls))

form.2.fs_gender <- formula(paste(endo.1, "~", inst.1, "+", time.controls, "+", demo.controls_gender))
form.3.fs_gender <- formula(paste(endo.1, "~", inst.1, "+", time.controls, "+", demo.controls_gender, "+", case.controls))

form.2_gender <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls_gender, "|", inst.1, "+", time.controls, "+", demo.controls_gender))
form.3_gender <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls_gender, "+", case.controls, "|", inst.1, "+", time.controls, "+", demo.controls_gender, "+", case.controls))

form.3.fs_drug <- formula(paste(endo.1, "~", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls_drug))
form.3_drug <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "+", case.controls_drug, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls_drug))

form.3.fs_prop <- formula(paste(endo.1, "~", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls_prop))
form.3_prop <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "+", case.controls_prop, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls_prop))

form.3.fs_weap <- formula(paste(endo.1, "~", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls_weapon))
form.3_weap <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "+", case.controls_weapon, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls_weapon))

form.1.fs_v <- formula(paste(endo.1, "~", inst.1, "+", time.controls_violent))
form.2.fs_v <- formula(paste(endo.1, "~", inst.1, "+", time.controls_violent, "+", demo.controls))
form.3.fs_v <- formula(paste(endo.1, "~", inst.1, "+", time.controls_violent, "+", demo.controls, "+", case.controls))

form.1_v <- formula(paste(outc.1, "~", endo.1, "+" , time.controls_violent, "|", inst.1, "+", time.controls_violent))
form.2_v <- formula(paste(outc.1, "~", endo.1, "+" , time.controls_violent, "+", demo.controls, "|", inst.1, "+", time.controls_violent, "+", demo.controls))
form.3_v <- formula(paste(outc.1, "~", endo.1, "+" , time.controls_violent, "+", demo.controls, "+", case.controls, "|", inst.1, "+", time.controls_violent, "+", demo.controls, "+", case.controls))

form.3.fs_prior <- formula(paste(endo.1, "~", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls_prior))
form.3_prior <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "+", case.controls_prior, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls_prior))

createTable <- function(split, cat,  m1, m2, m3, se, f1, f2, f3, n1, n2, n3){
  last.cases$split <- last.cases[[split]]
  c1 <- round(coeftest.cluster(last.cases[split == cat], m1, cluster1 = "judge_cat"),2)
  c1 <- c1[grep("judgeiv_hd", rownames(c1)), ]
  
  c2 <- round(coeftest.cluster(last.cases[split == cat], m2, cluster1 = "judge_cat"),2)
  c2 <- c2[grep("judgeiv_hd", rownames(c2)), ]
  
  c3 <- round(coeftest.cluster(last.cases[split == cat], m3, cluster1 = "judge_cat"),2)
  c3 <- c3[grep("judgeiv_hd", rownames(c3)), ]
  
  p1 <- data.table(cbind(c1["Estimate"], c2["Estimate"], c3["Estimate"]))
  p1$V4 <- 1
  p1$V5 <- 1:nrow(p1)
  p1$V6 <- "Estimate"
  
  p2 <- data.table(t(as.matrix(se)))
  colnames(p2) <- c("V1", "V2", "V3")
  p2$V4 <- 2
  p2$V5 <- 1:nrow(p2)
  p2$V6 <- "Std. Error"
  
  p3 <- data.table(round(cbind(f1, f2, f3), 2))
  colnames(p3) <- c("V1", "V2", "V3")
  p3$V4 <- 3
  p3$V5 <- 1:nrow(p3)
  p3$V6 <- "F-stat"
  
  table <- rbind(p1, p2, p3)
  table <- table[order(V5, V4), ]
  table$V4 <- table$V5 <- NULL
  p4 <- data.table(n1$nobs, n2$nobs, n3$nobs)    
  p4$V6 <- "N"
  
  table <- rbind(table, p4)
  table$V7 <- paste0(split,": ", cat)
  table <- table[, c("V7", "V6", "V1", "V2", "V3")]
  colnames(table) <- c("Subset", "Variable", "Model 1", "Model 2", "Model 3")
  
  return(table)
}

# -------------------
# Race: Black
# -------------------
m1d1.b <- lm(form.1.fs, data = last.cases[race4 == "B"])
m1d2.b <- lm(form.2.fs_race, data = last.cases[race4 == "B"])
m1d3.b <- lm(form.3.fs_race, data = last.cases[race4 == "B"])

m1a1.b <- ivreg(form.1, data = last.cases[race4 == "B"])
m1a2.b <- ivreg(form.2_race, data = last.cases[race4 == "B"])
m1a3.b <- ivreg(form.3_race, data = last.cases[race4 == "B"])

results.se1 <- list()
for(i in 1:500) {
  
  last.casesS <- last.cases[race4 == "B", ]
  out <- last.casesS[, this_id_fa[sample.int(.N, .N, TRUE)], by = c("judge_cat", "court_time1")]
  last.cases.boot <- last.casesS[last.casesS$this_id_fa %in% out$V1, ]

  m1d1.b0 <- lm(form.1.fs, data = last.cases.boot)
  m1d2.b0 <- lm(form.2.fs_race, data = last.cases.boot)
  m1d3.b0 <- lm(form.3.fs_race, data = last.cases.boot)
  
  results.se1[[i]] <- c(coefficients(m1d1.b0)["judgeiv_hd"],
                       coefficients(m1d2.b0)["judgeiv_hd"],
                       coefficients(m1d3.b0)["judgeiv_hd"]) 
}

results.se.f1 <- apply(do.call('rbind', results.se1), 2, sd)

m1a1_d.b <- summary(m1a1.b, vcov = sandwich, diagnostic = T)
m1a2_d.b <- summary(m1a2.b, vcov = sandwich, diagnostic = T)
m1a3_d.b <- summary(m1a3.b, vcov = sandwich, diagnostic = T)

m1a1_d1.b <- m1a1_d.b$diagnostics[grep("Weak", rownames(m1a1_d.b$diagnostics)), 3]
m1a2_d1.b <- m1a2_d.b$diagnostics[grep("Weak", rownames(m1a2_d.b$diagnostics)), 3]
m1a3_d1.b <- m1a3_d.b$diagnostics[grep("Weak", rownames(m1a3_d.b$diagnostics)), 3]

table.black <- createTable(split = "race4",cat = "B", 
                           m1 = m1d1.b, m2 = m1d2.b, m3 = m1d3.b, 
                           se = results.se.f1,
                           f1 = m1a1_d1.b, f2 = m1a2_d1.b, f3 = m1a3_d1.b,
                           n1 = m1a1.b, n2 = m1a2.b, n3 = m1a3.b)

# -------------------
# Race: White
# -------------------
m1d1.w <- lm(form.1.fs, data = last.cases[race4 == "W"])
m1d2.w <- lm(form.2.fs_race, data = last.cases[race4 == "W"])
m1d3.w <- lm(form.3.fs_race, data = last.cases[race4 == "W"])

results.se2 <- list()
for(i in 1:500) {
  
  last.casesS <- last.cases[race4 == "W", ]
  out <- last.casesS[, this_id_fa[sample.int(.N, .N, TRUE)], by = c("judge_cat", "court_time1")]
  last.cases.boot <- last.casesS[last.casesS$this_id_fa %in% out$V1, ]

  m1d1.w0 <- lm(form.1.fs, data = last.cases.boot)
  m1d2.w0 <- lm(form.2.fs_race, data = last.cases.boot)
  m1d3.w0 <- lm(form.3.fs_race, data = last.cases.boot)
  
  results.se2[[i]] <- c(coefficients(m1d1.w0)["judgeiv_hd"],
                        coefficients(m1d2.w0)["judgeiv_hd"],
                        coefficients(m1d3.w0)["judgeiv_hd"]) 
}

results.se.f2 <- apply(do.call('rbind', results.se2), 2, sd)

m1a1.w <- ivreg(form.1, data = last.cases[race4 == "W"])
m1a2.w <- ivreg(form.2_race, data = last.cases[race4 == "W"])
m1a3.w <- ivreg(form.3_race, data = last.cases[race4 == "W"])

m1a1_d.w <- summary(m1a1.w, vcov = sandwich, diagnostic = T)
m1a2_d.w <- summary(m1a2.w, vcov = sandwich, diagnostic = T)
m1a3_d.w <- summary(m1a3.w, vcov = sandwich, diagnostic = T)

m1a1_d1.w <- m1a1_d.w$diagnostics[grep("Weak", rownames(m1a1_d.w$diagnostics)), 3]
m1a2_d1.w <- m1a2_d.w$diagnostics[grep("Weak", rownames(m1a2_d.w$diagnostics)), 3]
m1a3_d1.w <- m1a3_d.w$diagnostics[grep("Weak", rownames(m1a3_d.w$diagnostics)), 3]

table.white <- createTable(split = "race4",cat = "W", 
                           m1 = m1d1.w, m2 = m1d2.w, m3 = m1d3.w, 
                           se = results.se.f2, 
                           f1 = m1a1_d1.w, f2 = m1a2_d1.w, f3 = m1a3_d1.w,
                           n1 = m1a1.w, n2 = m1a2.w, n3 = m1a3.w)

# -------------------
# Race: Hispanic
# -------------------
m1d1.h <- lm(form.1.fs, data = last.cases[race4 == "H"])
m1d2.h <- lm(form.2.fs_race, data = last.cases[race4 == "H"])
m1d3.h <- lm(form.3.fs_race, data = last.cases[race4 == "H"])

results.se3 <- list()
for(i in 1:500) {
  
  last.casesS <- last.cases[race4 == "H", ]
  out <- last.casesS[, this_id_fa[sample.int(.N, .N, TRUE)], by = c("judge_cat", "court_time1")]
  last.cases.boot <- last.casesS[last.casesS$this_id_fa %in% out$V1, ]

  m1d1.h0 <- lm(form.1.fs, data = last.cases.boot)
  m1d2.h0 <- lm(form.2.fs_race, data = last.cases.boot)
  m1d3.h0 <- lm(form.3.fs_race, data = last.cases.boot)
  
  results.se3[[i]] <- c(coefficients(m1d1.h0)["judgeiv_hd"],
                        coefficients(m1d2.h0)["judgeiv_hd"],
                        coefficients(m1d3.h0)["judgeiv_hd"]) 
}

results.se.f3 <- apply(do.call('rbind', results.se3), 2, sd)

m1a1.h <- ivreg(form.1, data = last.cases[race4 == "H"])
m1a2.h <- ivreg(form.2_race, data = last.cases[race4 == "H"])
m1a3.h <- ivreg(form.3_race, data = last.cases[race4 == "H"])

m1a1_d.h <- summary(m1a1.h, vcov = sandwich, diagnostic = T)
m1a2_d.h <- summary(m1a2.h, vcov = sandwich, diagnostic = T)
m1a3_d.h <- summary(m1a3.h, vcov = sandwich, diagnostic = T)

m1a1_d1.h <- m1a1_d.h$diagnostics[grep("Weak", rownames(m1a1_d.h$diagnostics)), 3]
m1a2_d1.h <- m1a2_d.h$diagnostics[grep("Weak", rownames(m1a2_d.h$diagnostics)), 3]
m1a3_d1.h <- m1a3_d.h$diagnostics[grep("Weak", rownames(m1a3_d.h$diagnostics)), 3]

table.hispanic <- createTable(split = "race4",cat = "H", 
                              m1 = m1d1.h, m2 = m1d2.h, m3 = m1d3.h,
                              se = results.se.f3,
                              f1 = m1a1_d1.h, f2 = m1a2_d1.h, f3 = m1a3_d1.h,
                              n1 = m1a1.h, n2 = m1a2.h, n3 = m1a3.h)
# -------------------
# Prior case: 1
# -------------------
m1d1.p <- lm(form.1.fs, data = last.cases[any_prior_case == 1])
m1d2.p <- lm(form.2.fs, data = last.cases[any_prior_case == 1])
m1d3.p <- lm(form.3.fs_prior, data = last.cases[any_prior_case == 1])

results.se4 <- list()
for(i in 1:500) {
  last.casesS <- last.cases[any_prior_case == 1, ]
  out <- last.casesS[, this_id_fa[sample.int(.N, .N, TRUE)], by = c("judge_cat", "court_time1")]
  last.cases.boot <- last.casesS[last.casesS$this_id_fa %in% out$V1, ]

  m1d1.p0 <- lm(form.1.fs, data = last.cases.boot)
  m1d2.p0 <- lm(form.2.fs, data = last.cases.boot)
  m1d3.p0 <- lm(form.3.fs_prior, data = last.cases.boot)
  
  results.se4[[i]] <- c(coefficients(m1d1.p0)["judgeiv_hd"],
                        coefficients(m1d2.p0)["judgeiv_hd"],
                        coefficients(m1d3.p0)["judgeiv_hd"]) 
}

results.se.f4 <- apply(do.call('rbind', results.se4), 2, sd)

m1a1.p <- ivreg(form.1, data = last.cases[any_prior_case == 1])
m1a2.p <- ivreg(form.2, data = last.cases[any_prior_case == 1])
m1a3.p <- ivreg(form.3_prior, data = last.cases[any_prior_case == 1])

m1a1_d.p <- summary(m1a1.p, vcov = sandwich, diagnostic = T)
m1a2_d.p <- summary(m1a2.p, vcov = sandwich, diagnostic = T)
m1a3_d.p <- summary(m1a3.p, vcov = sandwich, diagnostic = T)

m1a1_d1.p <- m1a1_d.p$diagnostics[grep("Weak", rownames(m1a1_d.p$diagnostics)), 3]
m1a2_d1.p <- m1a2_d.p$diagnostics[grep("Weak", rownames(m1a2_d.p$diagnostics)), 3]
m1a3_d1.p <- m1a3_d.p$diagnostics[grep("Weak", rownames(m1a3_d.p$diagnostics)), 3]

table.prior <- createTable(split = "any_prior_case",cat = 1, 
                           m1 = m1d1.p, m2 = m1d2.p, m3 = m1d3.p,
                           se = results.se.f4,
                           f1 = m1a1_d1.p, f2 = m1a2_d1.p, f3 = m1a3_d1.p,
                           n1 = m1a1.p, n2 = m1a2.p, n3 = m1a3.p)

# -------------------
# Prior case: 0
# -------------------
m1d1.np <- lm(form.1.fs, data = last.cases[any_prior_case == 0])
m1d2.np <- lm(form.2.fs, data = last.cases[any_prior_case == 0])
m1d3.np <- lm(form.3.fs_prior, data = last.cases[any_prior_case == 0])

results.se5 <- list()
for(i in 1:500) {
  last.casesS <- last.cases[any_prior_case == 0, ]
  out <- last.casesS[, this_id_fa[sample.int(.N, .N, TRUE)], by = c("judge_cat", "court_time1")]
  last.cases.boot <- last.casesS[last.casesS$this_id_fa %in% out$V1, ]

  m1d1.np0 <- lm(form.1.fs, data = last.cases.boot)
  m1d2.np0 <- lm(form.2.fs, data = last.cases.boot)
  m1d3.np0 <- lm(form.3.fs_prior, data = last.cases.boot)
  
  results.se5[[i]] <- c(coefficients(m1d1.np0)["judgeiv_hd"],
                        coefficients(m1d2.np0)["judgeiv_hd"],
                        coefficients(m1d3.np0)["judgeiv_hd"]) 
}

results.se.f5 <- apply(do.call('rbind', results.se5), 2, sd)

m1a1.np <- ivreg(form.1, data = last.cases[any_prior_case == 0])
m1a2.np <- ivreg(form.2, data = last.cases[any_prior_case == 0])
m1a3.np <- ivreg(form.3_prior, data = last.cases[any_prior_case == 0])

m1a1_d.np <- summary(m1a1.np, vcov = sandwich, diagnostic = T)
m1a2_d.np <- summary(m1a2.np, vcov = sandwich, diagnostic = T)
m1a3_d.np <- summary(m1a3.np, vcov = sandwich, diagnostic = T)

m1a1_d1.np <- m1a1_d.np$diagnostics[grep("Weak", rownames(m1a1_d.np$diagnostics)), 3]
m1a2_d1.np <- m1a2_d.np$diagnostics[grep("Weak", rownames(m1a2_d.np$diagnostics)), 3]
m1a3_d1.np <- m1a3_d.np$diagnostics[grep("Weak", rownames(m1a3_d.np$diagnostics)), 3]

table.noprior <- createTable(split = "any_prior_case",cat = 0, 
                             m1 = m1d1.np, m2 = m1d2.np, m3 = m1d3.np, 
                             se = results.se.f5,
                             f1 = m1a1_d1.np, f2 = m1a2_d1.np, f3 = m1a3_d1.np,
                             n1 = m1a1.np, n2 = m1a2.np, n3 = m1a3.np)
# -------------------
# Gender: M
# -------------------
m1d1.m <- lm(form.1.fs, data = last.cases[female == 0])
m1d2.m <- lm(form.2.fs_gender, data = last.cases[female == 0])
m1d3.m <- lm(form.3.fs_gender, data = last.cases[female == 0])

results.se6 <- list()
for(i in 1:500) {
  last.casesS <- last.cases[female == 0, ]
  out <- last.casesS[, this_id_fa[sample.int(.N, .N, TRUE)], by = c("judge_cat", "court_time1")]
  last.cases.boot <- last.casesS[last.casesS$this_id_fa %in% out$V1, ]
  
  m1d1.m0 <- lm(form.1.fs, data = last.cases.boot)
  m1d2.m0 <- lm(form.2.fs_gender, data = last.cases.boot)
  m1d3.m0 <- lm(form.3.fs_gender, data = last.cases.boot)
  
  results.se6[[i]] <- c(coefficients(m1d1.m0)["judgeiv_hd"],
                        coefficients(m1d2.m0)["judgeiv_hd"],
                        coefficients(m1d3.m0)["judgeiv_hd"]) 
}

results.se.f6 <- apply(do.call('rbind', results.se6), 2, sd)

m1a1.m <- ivreg(form.1, data = last.cases[female == 0])
m1a2.m <- ivreg(form.2_gender, data = last.cases[female == 0])
m1a3.m <- ivreg(form.3_gender, data = last.cases[female == 0])

m1a1_d.m <- summary(m1a1.m, vcov = sandwich, diagnostic = T)
m1a2_d.m <- summary(m1a2.m, vcov = sandwich, diagnostic = T)
m1a3_d.m <- summary(m1a3.m, vcov = sandwich, diagnostic = T)

m1a1_d1.m <- m1a1_d.m$diagnostics[grep("Weak", rownames(m1a1_d.m$diagnostics)), 3]
m1a2_d1.m <- m1a2_d.m$diagnostics[grep("Weak", rownames(m1a2_d.m$diagnostics)), 3]
m1a3_d1.m <- m1a3_d.m$diagnostics[grep("Weak", rownames(m1a3_d.m$diagnostics)), 3]

table.male <- createTable(split = "female",cat = 0, 
                          m1 = m1d1.m, m2 = m1d2.m, m3 = m1d3.m, 
                          se = results.se.f6,
                          f1 = m1a1_d1.m, f2 = m1a2_d1.m, f3 = m1a3_d1.m,
                          n1 = m1a1.m, n2 = m1a2.m, n3 = m1a3.m)

# -------------------
# Gender: F
# -------------------
m1d1.f <- lm(form.1.fs, data = last.cases[female == 1])
m1d2.f <- lm(form.2.fs_gender, data = last.cases[female == 1])
m1d3.f <- lm(form.3.fs_gender, data = last.cases[female == 1])

results.se7 <- list()
for(i in 1:500) {
  last.casesS <- last.cases[female == 1, ]
  out <- last.casesS[, this_id_fa[sample.int(.N, .N, TRUE)], by = c("judge_cat", "court_time1")]
  last.cases.boot <- last.casesS[last.casesS$this_id_fa %in% out$V1, ]
  
  m1d1.f0 <- lm(form.1.fs, data = last.cases.boot)
  m1d2.f0 <- lm(form.2.fs_gender, data = last.cases.boot)
  m1d3.f0 <- lm(form.3.fs_gender, data = last.cases.boot)
  
  results.se7[[i]] <- c(coefficients(m1d1.f0)["judgeiv_hd"],
                        coefficients(m1d2.f0)["judgeiv_hd"],
                        coefficients(m1d3.f0)["judgeiv_hd"]) 
}

results.se.f7 <- apply(do.call('rbind', results.se7), 2, sd)

m1a1.f <- ivreg(form.1, data = last.cases[female == 1])
m1a2.f <- ivreg(form.2_gender, data = last.cases[female == 1])
m1a3.f <- ivreg(form.3_gender, data = last.cases[female == 1])

m1a1_d.f <- summary(m1a1.f, vcov = sandwich, diagnostic = T)
m1a2_d.f <- summary(m1a2.f, vcov = sandwich, diagnostic = T)
m1a3_d.f <- summary(m1a3.f, vcov = sandwich, diagnostic = T)

m1a1_d1.f <- m1a1_d.f$diagnostics[grep("Weak", rownames(m1a1_d.f$diagnostics)), 3]
m1a2_d1.f <- m1a2_d.f$diagnostics[grep("Weak", rownames(m1a2_d.f$diagnostics)), 3]
m1a3_d1.f <- m1a3_d.f$diagnostics[grep("Weak", rownames(m1a3_d.f$diagnostics)), 3]

table.female <- createTable(split = "female",cat = 1, 
                          m1 = m1d1.f, m2 = m1d2.f, m3 = m1d3.f, 
                          se = results.se.f7,
                          f1 = m1a1_d1.f, f2 = m1a2_d1.f, f3 = m1a3_d1.f,
                          n1 = m1a1.f, n2 = m1a2.f, n3 = m1a3.f)

# -------------------
# Offense type: Drug
# -------------------
m1d1.d <- lm(form.1.fs, data = last.cases[any_drug == 1])
m1d2.d <- lm(form.2.fs, data = last.cases[any_drug == 1])
m1d3.d <- lm(form.3.fs_drug, data = last.cases[any_drug == 1])

results.se8 <- list()
for(i in 1:500) {
  last.casesS <- last.cases[any_drug == 1, ]
  out <- last.casesS[, this_id_fa[sample.int(.N, .N, TRUE)], by = c("judge_cat", "court_time1")]
  last.cases.boot <- last.casesS[last.casesS$this_id_fa %in% out$V1, ]
  
  m1d1.d0 <- lm(form.1.fs, data = last.cases.boot)
  m1d2.d0 <- lm(form.2.fs, data = last.cases.boot)
  m1d3.d0 <- lm(form.3.fs_drug, data = last.cases.boot)
  
  results.se8[[i]] <- c(coefficients(m1d1.d0)["judgeiv_hd"],
                        coefficients(m1d2.d0)["judgeiv_hd"],
                        coefficients(m1d3.d0)["judgeiv_hd"]) 
}

results.se.f8 <- apply(do.call('rbind', results.se8), 2, sd)

m1a1.d <- ivreg(form.1, data = last.cases[any_drug == 1])
m1a2.d <- ivreg(form.2, data = last.cases[any_drug == 1])
m1a3.d <- ivreg(form.3_drug, data = last.cases[any_drug == 1])

m1a1_d.d <- summary(m1a1.d, vcov = sandwich, diagnostic = T)
m1a2_d.d <- summary(m1a2.d, vcov = sandwich, diagnostic = T)
m1a3_d.d <- summary(m1a3.d, vcov = sandwich, diagnostic = T)

m1a1_d1.d <- m1a1_d.d$diagnostics[grep("Weak", rownames(m1a1_d.d$diagnostics)), 3]
m1a2_d1.d <- m1a2_d.d$diagnostics[grep("Weak", rownames(m1a2_d.d$diagnostics)), 3]
m1a3_d1.d <- m1a3_d.d$diagnostics[grep("Weak", rownames(m1a3_d.d$diagnostics)), 3]

table.drug <- createTable(split = "any_drug",cat = 1, 
                          m1 = m1d1.d, m2 = m1d2.d, m3 = m1d3.d, 
                          se = results.se.f8, 
                          f1 = m1a1_d1.d, f2 = m1a2_d1.d, f3 = m1a3_d1.d,
                          n1 = m1a1.d, n2 = m1a2.d, n3 = m1a3.d)
# -------------------
# Offense type: No Drug
# -------------------
m1d1.nd <- lm(form.1.fs, data = last.cases[any_drug == 0])
m1d2.nd <- lm(form.2.fs, data = last.cases[any_drug == 0])
m1d3.nd <- lm(form.3.fs_drug, data = last.cases[any_drug == 0])

results.se9 <- list()
for(i in 1:500) {
  last.casesS <- last.cases[any_drug == 0, ]
  out <- last.casesS[, this_id_fa[sample.int(.N, .N, TRUE)], by = c("judge_cat", "court_time1")]
  last.cases.boot <- last.casesS[last.casesS$this_id_fa %in% out$V1, ]
  
  m1d1.nd0 <- lm(form.1.fs, data = last.cases.boot)
  m1d2.nd0 <- lm(form.2.fs, data = last.cases.boot)
  m1d3.nd0 <- lm(form.3.fs_drug, data = last.cases.boot)
  
  results.se9[[i]] <- c(coefficients(m1d1.nd0)["judgeiv_hd"],
                        coefficients(m1d2.nd0)["judgeiv_hd"],
                        coefficients(m1d3.nd0)["judgeiv_hd"]) 
}

results.se.f9 <- apply(do.call('rbind', results.se9), 2, sd)

m1a1.nd <- ivreg(form.1, data = last.cases[any_drug == 0])
m1a2.nd <- ivreg(form.2, data = last.cases[any_drug == 0])
m1a3.nd <- ivreg(form.3_drug, data = last.cases[any_drug == 0])

m1a1_d.nd <- summary(m1a1.nd, vcov = sandwich, diagnostic = T)
m1a2_d.nd <- summary(m1a2.nd, vcov = sandwich, diagnostic = T)
m1a3_d.nd <- summary(m1a3.nd, vcov = sandwich, diagnostic = T)

m1a1_d1.nd <- m1a1_d.nd$diagnostics[grep("Weak", rownames(m1a1_d.nd$diagnostics)), 3]
m1a2_d1.nd <- m1a2_d.nd$diagnostics[grep("Weak", rownames(m1a2_d.nd$diagnostics)), 3]
m1a3_d1.nd <- m1a3_d.nd$diagnostics[grep("Weak", rownames(m1a3_d.nd$diagnostics)), 3]

table.ndrug <- createTable(split = "any_drug",cat = 0, 
                           m1 = m1d1.nd, m2 = m1d2.nd, m3 = m1d3.nd, 
                           se = results.se.f9, 
                           f1 = m1a1_d1.nd, f2 = m1a2_d1.nd, f3 = m1a3_d1.nd,
                           n1 = m1a1.nd, n2 = m1a2.nd, n3 = m1a3.nd)


# -------------------
# Offense type: Prop
# -------------------
m1d1.prop <- lm(form.1.fs, data = last.cases[any_prop == 1])
m1d2.prop <- lm(form.2.fs, data = last.cases[any_prop == 1])
m1d3.prop <- lm(form.3.fs_prop, data = last.cases[any_prop == 1])

results.se10 <- list()
for(i in 1:500) {
  last.casesS <- last.cases[any_prop == 1, ]
  out <- last.casesS[, this_id_fa[sample.int(.N, .N, TRUE)], by = c("judge_cat", "court_time1")]
  last.cases.boot <- last.casesS[last.casesS$this_id_fa %in% out$V1, ]
  
  m1d1.p0 <- lm(form.1.fs, data = last.cases.boot)
  m1d2.p0 <- lm(form.2.fs, data = last.cases.boot)
  m1d3.p0 <- lm(form.3.fs_prop, data = last.cases.boot)
  
  results.se10[[i]] <- c(coefficients(m1d1.p0)["judgeiv_hd"],
                        coefficients(m1d2.p0)["judgeiv_hd"],
                        coefficients(m1d3.p0)["judgeiv_hd"]) 
}

results.se.f10 <- apply(do.call('rbind', results.se10), 2, sd)

m1a1.prop <- ivreg(form.1, data = last.cases[any_prop == 1])
m1a2.prop <- ivreg(form.2, data = last.cases[any_prop == 1])
m1a3.prop <- ivreg(form.3_prop, data = last.cases[any_prop == 1])

m1a1_d.prop <- summary(m1a1.prop, vcov = sandwich, diagnostic = T)
m1a2_d.prop <- summary(m1a2.prop, vcov = sandwich, diagnostic = T)
m1a3_d.prop <- summary(m1a3.prop, vcov = sandwich, diagnostic = T)

m1a1_d1.prop <- m1a1_d.prop$diagnostics[grep("Weak", rownames(m1a1_d.prop$diagnostics)), 3]
m1a2_d1.prop <- m1a2_d.prop$diagnostics[grep("Weak", rownames(m1a2_d.prop$diagnostics)), 3]
m1a3_d1.prop <- m1a3_d.prop$diagnostics[grep("Weak", rownames(m1a3_d.prop$diagnostics)), 3]

table.prop <- createTable(split = "any_prop",cat = 1, 
                          m1 = m1d1.prop, m2 = m1d2.prop, m3 = m1d3.prop, 
                          se = results.se.f10, 
                          f1 = m1a1_d1.prop, f2 = m1a2_d1.prop, f3 = m1a3_d1.prop,
                          n1 = m1a1.prop, n2 = m1a2.prop, n3 = m1a3.prop)
# -------------------
# Offense type: No Prop
# -------------------
m1d1.nprop <- lm(form.1.fs, data = last.cases[any_prop == 0])
m1d2.nprop <- lm(form.2.fs, data = last.cases[any_prop == 0])
m1d3.nprop <- lm(form.3.fs_prop, data = last.cases[any_prop == 0])

results.se11 <- list()
for(i in 1:500) {
  last.casesS <- last.cases[any_prop == 0, ]
  out <- last.casesS[, this_id_fa[sample.int(.N, .N, TRUE)], by = c("judge_cat", "court_time1")]
  last.cases.boot <- last.casesS[last.casesS$this_id_fa %in% out$V1, ]
  
  m1d1.np0 <- lm(form.1.fs, data = last.cases.boot)
  m1d2.np0 <- lm(form.2.fs, data = last.cases.boot)
  m1d3.np0 <- lm(form.3.fs_prop, data = last.cases.boot)
  
  results.se11[[i]] <- c(coefficients(m1d1.np0)["judgeiv_hd"],
                         coefficients(m1d2.np0)["judgeiv_hd"],
                         coefficients(m1d3.np0)["judgeiv_hd"]) 
}

results.se.f11 <- apply(do.call('rbind', results.se11), 2, sd)

m1a1.nprop <- ivreg(form.1, data = last.cases[any_prop == 0])
m1a2.nprop <- ivreg(form.2, data = last.cases[any_prop == 0])
m1a3.nprop <- ivreg(form.3_prop, data = last.cases[any_prop == 0])

m1a1_d.nprop <- summary(m1a1.nprop, vcov = sandwich, diagnostic = T)
m1a2_d.nprop <- summary(m1a2.nprop, vcov = sandwich, diagnostic = T)
m1a3_d.nprop <- summary(m1a3.nprop, vcov = sandwich, diagnostic = T)

m1a1_d1.nprop <- m1a1_d.nprop$diagnostics[grep("Weak", rownames(m1a1_d.nprop$diagnostics)), 3]
m1a2_d1.nprop <- m1a2_d.nprop$diagnostics[grep("Weak", rownames(m1a2_d.nprop$diagnostics)), 3]
m1a3_d1.nprop <- m1a3_d.nprop$diagnostics[grep("Weak", rownames(m1a3_d.nprop$diagnostics)), 3]

table.nprop <- createTable(split = "any_prop",cat = 0, 
                           m1 = m1d1.nprop, m2 = m1d2.nprop, m3 = m1d3.nprop, 
                           se = results.se.f11, 
                           f1 = m1a1_d1.nprop, f2 = m1a2_d1.nprop, f3 = m1a3_d1.nprop,
                           n1 = m1a1.nprop, n2 = m1a2.nprop, n3 = m1a3.nprop)

# -------------------
# Offense type: Weap
# -------------------
m1d1.weap <- lm(form.1.fs, data = last.cases[any_weapon == 1])
m1d2.weap <- lm(form.2.fs, data = last.cases[any_weapon == 1])
m1d3.weap <- lm(form.3.fs_weap, data = last.cases[any_weapon == 1])

results.se12 <- list()
for(i in 1:500) {
  
  last.casesS <- last.cases[any_weapon == 1, ]
  out <- last.casesS[, this_id_fa[sample.int(.N, .N, TRUE)], by = c("judge_cat", "court_time1")]
  last.cases.boot <- last.casesS[last.casesS$this_id_fa %in% out$V1, ]
  
  m1d1.w0 <- lm(form.1.fs, data = last.cases.boot)
  m1d2.w0 <- lm(form.2.fs, data = last.cases.boot)
  m1d3.w0 <- lm(form.3.fs_weap, data = last.cases.boot)
  
  results.se12[[i]] <- c(coefficients(m1d1.w0)["judgeiv_hd"],
                         coefficients(m1d2.w0)["judgeiv_hd"],
                         coefficients(m1d3.w0)["judgeiv_hd"]) 
}

results.se.f12 <- apply(do.call('rbind', results.se12), 2, sd)

m1a1.weap <- ivreg(form.1, data = last.cases[any_weapon == 1])
m1a2.weap <- ivreg(form.2, data = last.cases[any_weapon == 1])
m1a3.weap <- ivreg(form.3_weap, data = last.cases[any_weapon == 1])

m1a1_d.weap <- summary(m1a1.weap, vcov = sandwich, diagnostic = T)
m1a2_d.weap <- summary(m1a2.weap, vcov = sandwich, diagnostic = T)
m1a3_d.weap <- summary(m1a3.weap, vcov = sandwich, diagnostic = T)

m1a1_d1.weap <- m1a1_d.weap$diagnostics[grep("Weak", rownames(m1a1_d.weap$diagnostics)), 3]
m1a2_d1.weap <- m1a2_d.weap$diagnostics[grep("Weak", rownames(m1a2_d.weap$diagnostics)), 3]
m1a3_d1.weap <- m1a3_d.weap$diagnostics[grep("Weak", rownames(m1a3_d.weap$diagnostics)), 3]

table.weap <- createTable(split = "any_weapon",cat = 1, 
                          m1 = m1d1.weap, m2 = m1d2.weap, m3 = m1d3.weap, 
                          se = results.se.f12, 
                          f1 = m1a1_d1.weap, f2 = m1a2_d1.weap, f3 = m1a3_d1.weap,
                          n1 = m1a1.weap, n2 = m1a2.weap, n3 = m1a3.weap)

# -------------------
# Offense type: No Weap
# -------------------
m1d1.nweap <- lm(form.1.fs, data = last.cases[any_weapon == 0])
m1d2.nweap <- lm(form.2.fs, data = last.cases[any_weapon == 0])
m1d3.nweap <- lm(form.3.fs_weap, data = last.cases[any_weapon == 0])

results.se13 <- list()
for(i in 1:500) {
  
  last.casesS <- last.cases[any_weapon == 0, ]
  out <- last.casesS[, this_id_fa[sample.int(.N, .N, TRUE)], by = c("judge_cat", "court_time1")]
  last.cases.boot <- last.casesS[last.casesS$this_id_fa %in% out$V1, ]
  
  m1d1.nw0 <- lm(form.1.fs, data = last.cases.boot)
  m1d2.nw0 <- lm(form.2.fs, data = last.cases.boot)
  m1d3.nw0 <- lm(form.3.fs_weap, data = last.cases.boot)
  
  results.se13[[i]] <- c(coefficients(m1d1.nw0)["judgeiv_hd"],
                         coefficients(m1d2.nw0)["judgeiv_hd"],
                         coefficients(m1d3.nw0)["judgeiv_hd"]) 
}

results.se.f13 <- apply(do.call('rbind', results.se13), 2, sd)

m1a1.nweap <- ivreg(form.1, data = last.cases[any_weapon == 0])
m1a2.nweap <- ivreg(form.2, data = last.cases[any_weapon == 0])
m1a3.nweap <- ivreg(form.3_weap, data = last.cases[any_weapon == 0])

m1a1_d.nweap <- summary(m1a1.nweap, vcov = sandwich, diagnostic = T)
m1a2_d.nweap <- summary(m1a2.nweap, vcov = sandwich, diagnostic = T)
m1a3_d.nweap <- summary(m1a3.nweap, vcov = sandwich, diagnostic = T)

m1a1_d1.nweap <- m1a1_d.nweap$diagnostics[grep("Weak", rownames(m1a1_d.nweap$diagnostics)), 3]
m1a2_d1.nweap <- m1a2_d.nweap$diagnostics[grep("Weak", rownames(m1a2_d.nweap$diagnostics)), 3]
m1a3_d1.nweap <- m1a3_d.nweap$diagnostics[grep("Weak", rownames(m1a3_d.nweap$diagnostics)), 3]

table.nweap <- createTable(split = "any_weapon",cat = 0, 
                           m1 = m1d1.nweap, m2 = m1d2.nweap, m3 = m1d3.nweap, 
                           se = results.se.f13, 
                           f1 = m1a1_d1.nweap, f2 = m1a2_d1.nweap, f3 = m1a3_d1.nweap,
                           n1 = m1a1.nweap, n2 = m1a2.nweap, n3 = m1a3.nweap)

# -------------------
# Offense type: Violent
# -------------------
m1d1.v <- lm(form.1.fs_v, data = last.cases[any_violence == 1])
m1d2.v <- lm(form.2.fs_v, data = last.cases[any_violence == 1])
m1d3.v <- lm(form.3.fs_v, data = last.cases[any_violence == 1])

results.se14 <- list()
for(i in 1:500) {
  
  last.casesS <- last.cases[any_violence == 1, ]
  out <- last.casesS[, this_id_fa[sample.int(.N, .N, TRUE)], by = c("judge_cat", "court_time1")]
  last.cases.boot <- last.casesS[last.casesS$this_id_fa %in% out$V1, ]
  
  m1d1.v0 <- lm(form.1.fs_v, data = last.cases.boot)
  m1d2.v0 <- lm(form.2.fs_v, data = last.cases.boot)
  m1d3.v0 <- lm(form.3.fs_v, data = last.cases.boot)
  
  results.se14[[i]] <- c(coefficients(m1d1.v0)["judgeiv_hd"],
                         coefficients(m1d2.v0)["judgeiv_hd"],
                         coefficients(m1d3.v0)["judgeiv_hd"]) 
}

results.se.f14 <- apply(do.call('rbind', results.se14), 2, sd)

m1a1.v <- ivreg(form.1_v, data = last.cases[any_violence == 1])
m1a2.v <- ivreg(form.2_v, data = last.cases[any_violence == 1])
m1a3.v <- ivreg(form.3_v, data = last.cases[any_violence == 1])

m1a1_d.v <- summary(m1a1.v, vcov = sandwich, diagnostic = T)
m1a2_d.v <- summary(m1a2.v, vcov = sandwich, diagnostic = T)
m1a3_d.v <- summary(m1a3.v, vcov = sandwich, diagnostic = T)

m1a1_d1.v <- m1a1_d.v$diagnostics[grep("Weak", rownames(m1a1_d.v$diagnostics)), 3]
m1a2_d1.v <- m1a2_d.v$diagnostics[grep("Weak", rownames(m1a2_d.v$diagnostics)), 3]
m1a3_d1.v <- m1a3_d.v$diagnostics[grep("Weak", rownames(m1a3_d.v$diagnostics)), 3]

table.v <- createTable(split = "any_violence",cat = 1, 
                       m1 = m1d1.v, m2 = m1d2.v, m3 = m1d3.v, 
                       se = results.se.f14, 
                       f1 = m1a1_d1.v, f2 = m1a2_d1.v, f3 = m1a3_d1.v,
                       n1 = m1a1.v, n2 = m1a2.v, n3 = m1a3.v)

# -------------------
# Offense type: Non-Violent
# -------------------
m1d1.nv <- lm(form.1.fs_v, data = last.cases[any_violence == 0])
m1d2.nv <- lm(form.2.fs_v, data = last.cases[any_violence == 0])
m1d3.nv <- lm(form.3.fs_v, data = last.cases[any_violence == 0])

results.se15 <- list()
for(i in 1:500) {
  
  last.casesS <- last.cases[any_violence == 0, ]
  out <- last.casesS[, this_id_fa[sample.int(.N, .N, TRUE)], by = c("judge_cat", "court_time1")]
  last.cases.boot <- last.casesS[last.casesS$this_id_fa %in% out$V1, ]
  
  m1d1.nv0 <- lm(form.1.fs_v, data = last.cases.boot)
  m1d2.nv0 <- lm(form.2.fs_v, data = last.cases.boot)
  m1d3.nv0 <- lm(form.3.fs_v, data = last.cases.boot)
  
  results.se15[[i]] <- c(coefficients(m1d1.nv0)["judgeiv_hd"],
                         coefficients(m1d2.nv0)["judgeiv_hd"],
                         coefficients(m1d3.nv0)["judgeiv_hd"]) 
}

results.se.f15 <- apply(do.call('rbind', results.se15), 2, sd)

m1a1.nv <- ivreg(form.1_v, data = last.cases[any_violence == 0])
m1a2.nv <- ivreg(form.2_v, data = last.cases[any_violence == 0])
m1a3.nv <- ivreg(form.3_v, data = last.cases[any_violence == 0])

m1a1_d.nv <- summary(m1a1.nv, vcov = sandwich, diagnostic = T)
m1a2_d.nv <- summary(m1a2.nv, vcov = sandwich, diagnostic = T)
m1a3_d.nv <- summary(m1a3.nv, vcov = sandwich, diagnostic = T)

m1a1_d1.nv <- m1a1_d.nv$diagnostics[grep("Weak", rownames(m1a1_d.nv$diagnostics)), 3]
m1a2_d1.nv <- m1a2_d.nv$diagnostics[grep("Weak", rownames(m1a2_d.nv$diagnostics)), 3]
m1a3_d1.nv <- m1a3_d.nv$diagnostics[grep("Weak", rownames(m1a3_d.nv$diagnostics)), 3]

table.nv <- createTable(split = "any_violence",cat = 0, 
                        m1 = m1d1.nv, m2 = m1d2.nv, m3 = m1d3.nv, 
                        se = results.se.f15, 
                        f1 = m1a1_d1.nv, f2 = m1a2_d1.nv, f3 = m1a3_d1.nv,
                        n1 = m1a1.nv, n2 = m1a2.nv, n3 = m1a3.nv)


all <- rbind(table.black, table.white, table.hispanic, table.prior, table.noprior, 
             table.male,table.female, 
             table.drug,table.ndrug,
             table.prop, table.nprop,
             table.weap, table.nweap,
             table.v, table.nv)
all

